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We study the statistics of encounters of Levy flights by introducing the concept of vicious Levy 
flights - distinct groups of walkers performing independent Levy flights with the process terminating 
upon the first encounter between walkers of different groups. We show that the probability that 
the process survives up to time t decays as t~" at late times. We compute a up to the second 
order in e-expansion, where e = ct — d, cr is the Levy exponent and d is the spatial dimension. For 
d = a, we find the exponent of the logarithmic decay exactly. Theoretical values of the exponents 
I*^ . are confirmed by numerical simulations. 
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*vj ' Diffusive processes with long range jumps play an important role in many physical, chemical and biological phenom- 

ena. A Levy flight is an example of such a process where the probability distribution of the length of an individual 
step, r, is governed by the power-law r~'''~°', where d is the dimension of the space and a is the Levy exponent. 
Smaller values of a therefore produce longer range jumps while for cr > 2, the mean jump length is finite and simple 
^ ■ diffusive behavior is recovered. Levy flights have been used to describe a wide range of processes including epidemic 
D ' spreading, transcription factor proteins binding to DNA, kinetic Ising models with long-range interactions, foraging 
H ■ animals and light propagation in disordered optical materials [iHg. While individual Levy flights have been studied 
' ' in great detail, the same is not true if we consider several distinct groups of Levy flights. One could, for example, be 
C^ ■ interested in the statistics of encounters between members of different groups. This question is relevant for processes 
c/2 ] where the outcome depends on the occurrence of such encounters. Examples include sharks and other marine animals 
searching for prey [7|, chemical reactions in turbulent environments [8|, electron- hole recombination in disordered 
media [9| and even male spider-monkeys encountering their mates or other aggressive males in the forest [lO| . 
Ch In this Letter we compute the survival probability, i.e. the probability that no two members of different groups 

——k of Levy flights have met up to time t. For the case of simple diffusion with exactly one particle in each group, this 
r^ ■ corresponds to the classic problem of Gaussian vicious walks ll|, i.e. walks that are prohibited from being on the 
Q \ same site at the same time, but remain independent otherwise. Here we generalize this concept to groups of Levy 
, O, flights under the same constraints. We term them vicious Levy flights (VLF). We consider p sets of particles with 
rii particles in each set, i — 1 . . .p, that are driven by Levy noise on the d dimensional regular lattice. A pairwise 
7—i • interset short-range (delta-function) interaction is introduced to guarantee that trajectories which continue beyond 
^ ' an intersection are discarded, i.e. have zero statistical weight. This terminates the process at the first encounter 
^!~ ' between members of different groups. Particles belonging to the same set do not interact. We note that Levy flights 
^^ ' are allowed to jump over each other, unlike ordinary random walks which can only jump to neighboring sites and can 
rvq ' not intersect with the vicious constraint. In d = 1 this means that the ordering is preserved for vicious walks but 

J • I not for VLF. For simplicity we assume that Levy exponents for all flights are the same. Generalization to the case of 

f- >) • different Levy exponents will be done elsewhere. At time t = all particles start in the vicinity of the origin. We are 

(^ ' interested in the survival probability of this system at late times. 

^-H . We start with a field theoretic formulation of the problem. Methods to formulate field theories for such stochastic 



systems are well established [12|, [13j. Specifically for Gaussian vicious walks, such a formulation exists [IJ] and the 
form of the action is known. We can adapt the action to our case by replacing the Laplacian V^ with the operator 
V^ that generates long-range jumps. This gives 

r ^ 

'5(0»,</'I)= /did''a:^[(/.lat(^,;+0lV>,]+ ^ X,j4{t,x)Mt,x)c^]it,x)<f>j{t,x). (1) 

where (j>i{x,t) are p complex order parameters corresponding to p different sets of equivalent Levy flights and A^^ 
are coupling constants corresponding to interset interactions. The non-intersection property of VLF arises from the 
choice Xij — > oo but we will show that to leading order the survival probability does not depend on the particular 
value of these coupling constants. This action is also similar to the action for the reaction-diffusion problem with 
long-range interactions 15|, Il6| . Power counting shows that the upper critical dimension for the above field theory is 



dc = <J for a < 2 and dc = 2 ior a > 2. VLF exhibit different phases (see inset Figure 1) depending on the values of 
d and a. In the mean field phase for d > dc (Region I), the survival probability of VLF is non-zero at infinite time 
because the walks become non-recurrent and particles can avoid each other for all time. For a > 2 (Region II) VLF 
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FIG. 1: a. In(S') vs ln(f) for two identical VLF in d = 1. a values from top tp bottom are 1.1, 1.3, 1.5, 1.7, 2.5 respectively. 
Symbols represent simulation data and solid lines correspond to best fit lines to the late time data. Inset: Domains of VLF 
exponents in the a — d plane. 

reproduce Gaussian vicious walks. For d < a < 2 (Region III) we expect fluctuations to play an important role. In 
this phase, we will obtain the critical behavior of the survival probability using e-expansion (e = a — d) around mean 
field theory in two-loop approximation. 

We now turn to the renormalization group analysis. The propagator given by ([T]) is Ta ' (s, k) ~ {s + k'^)^^. The 
particular form of the vertex in ([T]) leads to the fact that there are no diagrams which dress the propagator. This 
implies that the bare propagator is the exact propagator for the theory. The proper vertex is defined by factoring 
out external legs from the ordinary 4-point Green's function of ^ Gl,' '{si, kf, A). Here A = {\ij} is the collection of 
coupling constants and (sj, fc;) for I — 1 . . .A are the energy (Laplace image of time) and momenta respectively 17 1. 
This yields 



rlf'\s,,kr,\) = 



Gf/\si,ki;X) 

nl=ir(i'i)(s™,fc™) 



(2) 



The renormalized coupling, Xrij, is the value of the proper vertex at (s; — fi,ki — 0), for all I, where /i is a 
renormalization group flow scaling parameter. It is possible to sum all diagrams in the series with the result 



^Rij — Xij{l + Xijil) , 

where /i = {2n)^'^ J d'^q{2pL + ^q"^)^^. The value of the integral /i is given by /i = A/i^'/'^e"^, where 

eT{d/a)V{e/G) 2-''tt-''/'^ 
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is the geometric factor at the leading order in e — a -- d. By introducing the redefined coupling constant gjnj 
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we obtain the following renormalization group flow equations (see Appendix A for details): 

^~^ = (-£ + AgRij)gRij/a. 



(5) 



The fixed point of this flow is g* — eA^^. We note that this value of the fixed point is exact to all orders 
since all diagrams were taking into account in ([3]). The stability of the fixed point follows from the fact that 
—dl3ij/dgRij\gj^..=g/A = —e/a < in the VLF region, where Pij = ^dgmj/d^ is renormalization group beta function. 
We now consider the survival probability which is defined as the correlation function |14| 
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with the measure j2?0'''2?(/)exp[— 5]. The Feynman diagram of (jS]) at zero order is a vertex with 2N external legs. 
Similar to the case of the (2, 2)-vertex it is convenient to work with the truncated correlation function T{si, kf, A) = 
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FIG. 2: a. 1-loop diagram corresponding to the integral Ii b., c. 2-Ioop intgrals corresponding to /^ and I2 respectively. 

S{si, ki; ^)/Y[m=i ^ i^m, km). The finite renormalized truncated correlation function T]i{si, fc;; Aijj/i), wliere Xji = 
{^Rij} is a collection of renormalized coupling constants, is related to tlie bare truncated correlation function by 
Tf{{si, ki] Afl, fi) = Z(X, fj,)T{si, ki] A), where Z(\^ /i) is the scaling function. From this one obtains the renormalization 
group equation for T ii(si, kf, A/j, /j.) using the chain rule 

d d 

{^J■Tr- +/3»i^ h7)ri?(si,fc(;Afl,pi) = 0, (7) 

Cm cigRij 

where 7 = fidlnZ/dfi. At the fixed point ([7]) reduces to {d/d\n{fi) + -f^)r ]i{s , fi) — whose solution is 

Tnr^eMf%*diHp'))), (8) 

Jo 

where 7* = fidln Z/dfi\gj^-.^g,^. Since 7, is constant at the fixed point we have Tu ^ IJ-^'^* ■ The fact that the 
dimensions of field and the action are [t/i^] — [4>\ — k'^^^ and [S] ~ 1 implies that [5*] = 1. Thus it follows that the 
survival probability can only be a function of the dimensionless product [it. From this one infers that the asymptotic 
behavior of the survival probability is S'(i) ~ t~^* which gives a = 7*. In order to find Z one uses a normahzation 
condition on F^j that fixes the value of Z . This can be chosen as Ffl(s;, fcj; Ajj, [i) — \ when si = [i and fc; = for all 
I. This implies that Z = F(/^, 0;A)~^. F(/i,0;A) can be expressed as a series, up to two-loop order, of the integrals 
corresponding to the diagrams shown in Fig. 2 with appropriate combinatorial factors originating from the number 
of distinct ways in which propagators can be assigned to the same diagram, /i was evaluated before. The integral 
corresponding to the diagram 2c is 

^^J {2[i + 2k''){3fi + q'' + k'' + \k + q\'') ^' 

This can be evaluated using Mellin-Barnes representation [18] which replaces the sum in the denominator of |fc + q]"^ 
and g'^ by the product of these terms raised to some power. The result for I2 up to the leadiirg order reads 



,-^-/-rjL+^(-i°g(3/4)/4+c) (^^^ 



" r(o-/2)2 ' y2e^ ae 

where C = [ilj^-^^a/2) + log(47r)]/2 and ip'-'^^x) is standard digamma function. The details of the calculations are 
summarized in Appendix B. Knowing /i, I2 and the appropriate combinatorial factors allows us to evaluate F(/i, 0; A) 
and therefore Z{p,,X) as a series. Differentiating InZ with respect to /i and substituting A with A^, (inverting ([3])) 
and then taking the value at the fixed point gives us the survival probabiUty exponent a = 7*, with the final result 
(see Appendix C): 

a= ^ n^nj-e/cr + hl(3/4)Q£Vc^^ (H) 

l<i<j<p 

where Q — Q X]i<i<7<fe<i) ^^i'^j'^fc + X]i<i<i<T) ^i^ji''^i + ^^j ~ 2). At the critical dimension d — dc — u we see that the 
fixed point coincides with the Gaussian point. The interaction becomes marginal in the renormalization group sense. 
Equation ([S]) then yields the fiow equations for the running coupling constant 

xdgij {x)/dx = Agf^ {x)/a; (12) 

with initial condition gij{l) — gmj. Solving this and substituting the result into ([5]) we get S{t) — (Int)""', where 

ai = ^ Uirij. (13) 

l<i<j<p 
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FIG. 3: a as a function of a for d = 1. Symbols with error bars represent simulation data corresponding from top to bottom to 
A'^ = 2, 3, 4 VLF respectively. Lines correspond to 1-loop approximation from formula dill) for 1 < a < 2. For a < 1, a >2 lines 
represent the mean field and Gaussian exponents respectively. Inset: Same simulation data compared to 2-loop approximation 
from {ni) 

Now we describe the details of the numerical simulation that we used to confirm our results in d = 1. At i = we 
start with TV = X^iLi '^i particles belonging to p distinct sets placed equidistantly on the lattice. At each time step we 
generate N random variables, Xj, drawn from the uniform distribution on the interval (0, 1). Each particle jumps a 



distance U 



.-1/'^ 



with equal probability to the left or to the right. This procedure generates an independent Levy 
flight trajectory for each particle. The process stops whenever particles from different sets land on the same site. We 
perform ~ 10^ iterations for each set of parameter values. The survival probability S{t) is defined by the number 
of processes that survived beyond time t divided by the total number of iterations. Figure 1 shows the plot of the 
survival probability as a function of time for TV = 2 for different values of a. It is clear that at late times In S{t) is 
linear in Int verifying our predicted power-law decay S{t) ^ t^". The critical exponent a is evaluated from the slope 
of the best fit line of the late time data. 

We first consider systems with exactly one particle in each set. Figure 3 shows the value of exponent a for various 
values of a with the total number of particles A'' = 2, 3 and 4 in d = 1. Values of cr > 2 will reproduce Gaussian vicious 
walks and we therefore expect our exponents to approach the exact Fisher exponents llj as seen in Figure 3. For 
two VLF higher loop corrections are absent (see (fTTjl ) and the one- loop result is an exact result in agreement with the 
simulation. It is interesting to note that the survival probability for the N = 2 case is equivalent to the first return 
probability of a single Levy flight to the origin after time t which scales as ^-1+'^/'^ [ig| and matches our results. For 
(7 < 1, or d > dc, we expect mean fleld behavior where survival probability at late times approaches a non-zero value 
implying that there is a finite probability that Levy flights with a < 1 will never find each other. The fact that the 
survival probability tends to a constant at late times is reflected in the small values of a for a = 0.9. For I < a < 2, 
the mean fleld behavior is incorrect and we expect the fluctuations to shift the decay exponent to some non-trivial 
value. For a close to one, the 1-loop result is in good agreement with the simulation. For larger values of a, the 2-loop 
corrections perform better (see Fig. 3 inset). It is to be noted that the discrepancy between theory and simulation 
becomes large for higher values of A^ because the combinatorial factors in pT|) become large and we therefore need to 
keep higher order terms in e for the same degree of accuracy. It is interesting that the 1-loop approximation works 
reasonably well over the entire range of 1 < cr < 2 simply because the 1-loop term in (fTTj) happens to give the exact 
Fisher result if we set a = 2. We notice that in all cases the value of the survival probability exponent a increases 
with a starting from a ^ at a < 1 and rising to the value of the Fisher exponents for the equivalent Gaussian vicious 
walks. This is in contrast to diffusion-annihilation reactions with long-range jumps where the density of reactants 
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FIG. 4: a vs a for predator and prey problem in d = 1. Symbols represent simulation data for 4 predators and 1 prey (circles) 
and 3 predators and 2 prey (squares). Lines are 2-loop approximation from formula (|11[) . 

decays faster for smaller values of a jl6| . 

We now consider a system that consists of 2 sets of identical VLF with different numbers of independent particles in 
each set. We shall call one set predators and the other set prey. Figure 4 compares the values of the 2-loop exponents 
to the simulation results for various values of a for two different cases: 4 predators - 1 prey and 3 predators - 2 
prey. Similar to the previous case we have mean field and Gaussian behaviors for a < 1 and a > 2 respectively. The 
Gaussian case is also known as the lion-lamb problem and has been studied before [20!] ■ Unlike the lion-lamb problem, 
however, our results do not depend on the initial ordering of predators and prey because ordering is not preserved for 
VLF. For a given a and total number of predators plus prey, the number of potentially lethal encounters is maximized 
when the difference between the number of predators and prey is the smallest implying that the survival probability 
will decay faster as seen in Fig. 4. 

Our results suggest that it is interesting to solve the problem in the general case of particles with different diffusion 
constants and Levy exponents. The predictive power of the e-expansion for VLF, that we have demonstrated, should 
be useful in many applications of practical importance. Examples include the optimization of the predator-prey 
search 2l| or trapping probabilities 22|. Generalization to the case of intelligent predators, i.e. interacting with a 



prey by means of the long-range potential, may lead to different critical behavior 23l-l25| . Simple diffusion processes 
in power-law small world networks are effectively Levy flights [26] with the exponent a controlling the distribution 
of long-range links. Our work could be used to understand what network structure, or what a, would optimize the 
search and how much more efficient several independent searchers will be. 

The authors would like to acknowledge UC Merced start-up funds and a James S. McDonnell Foundation Award 
for Studying Complex Systems. 
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Here we derive the 1-loop integral 



Ji = /i(;i) = {2tt)-'^ / d''k(2fi + 2k'') 



(14) 



We will use dimensional regularization. First we notice that there is no angle dependence under the integral thus one 
can integrate out d — 1 angle variable and use alpha representation: 

(•+00 
X-^=r{X)-^ daa^-^ expi'-aX) (15) 

to handle Id momenta integral: 

+ CXJ +00 +00 +00 

/i = ^ / ^ = ^ / / dadyfcfe'^-iexp(-a/x-afc") = -^r(d/a) / daa"'^/" exp(-a/^), (16) 

2 J ^ + k<^ 2 J J 2(7 J 







where Sd = 2tt'^/'^ /T{d/2) is the area of the d-dimcnsional unit sphere. After taking the integral over a one has 

h = ^T{d|a)T{e|a)^i-''- = Ay.-^l'^e-^ + 0(e°), (17) 

la 

where A has been defined by the formula (|4]). 

Now we show details of deriving renormalization group flow equations ([5]). We start with equation ([3]) and express 
Aij in terms of A/j^j . The result reads 

^- = 1 ^T T ■ (18) 

t — Aflijii 

Multiplying left and right hand side of the last equation on pT^^" and redefining the coupling constant g^ij — 
ji^'^i'^Xiiij we infer that 

ix'^'^X,, = gmj/[l - gmjAe^^). (19) 

Differentiating left and right hand side of p^ with fi-M- we obtain 

(-e/a)/.-^/^A,, = -A,g^,V(g^i. - Ae-'f (20) 

Now we substitute (1191) into (j20l) and find beta function up to second order in small e and gmj expansion: 

P^3 = i-£9mj + Agl,^)/a + 0{eg^) (21) 

APPENDIX B 

Here we derive the 2-loop integral 

h ^ Hti) ^ {2^T^'^ f d'^kd'^q[{2fi + 2k^){3^i + k^+q^ + \k + qn]-\ (22) 



The term |fc + gj*^ leads to the appearance of angle integration. Nevertheless it is possible to avoid angle integration. 
The key idea is to use Mellin-Barnes representation [18|: 

1 _ /■+^- dz Y^ r(A + z)r(-z) ^23) 



(X + Y)^ J_,^ 2niX^+^ r(A) 

Applying MB formula twice we split the sum of two terms containing q integration into the factor of these terms 
raised to some power: 

+200 

d'^kd'^q f dz T{l + zi)T{-zi){3fi + k'' + q'^y^ 



2(27r)2'i J 2m fi + k'^ |fc + 9|'^(i+^i) 

— ioo 
-{-ioo 

d'^kd'^q f dzidz2T{i + zi)r{-zi + z2)r{^z2) {iti + k-^y^ 

(24) 



2(27r)2rf 7 {2mY ^i + k'' |fc + g|'^(i+^i)q'^(-^i+^2) ' 
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Now integral over q becomes standard: 

dS d/2, d-2(a,+a.) r(ai + 02 - d/2)V{d/2 - a,)T(d/2 ~ a^) 

(q2)ai((fc + g)2)a. r(ai)r(a2)r(d- 01-02) 

where ai — a{^zi + ^2)/2 and 02 = cr(l + zi)/2. Thus we wiU be left with integral over k of the form 



f d'^kk-^'''''^{i^jL + k'') 
''^ J 2n + 2k'' 



(26) 



The function under the integral does not depend on the angle and therefore Ik can be cast into one dimensional 
integral over momenta: 

4 = |i/V-"-ffii±fi:i (27) 

2a J fi + k 



We will compute this integral using alpha representation. 

+00 
Sri f a~^^~ k^ 1"^'^^ 

Ik = — / dkdaida2 -, r exp(— 3//ai — aik — Q;2/i — 0^2^) (28) 

2a- J ^{-Z2) 



After momenta integration we obtain 

+00 

h = ^'^^(l - 2£/^ - ^2) f ^^^^^^(^^ + a2f'/^+^'-'ar'-' exp(-3/xai - ^2^) (29) 

zal (-Z2) J 



First we will take care the integral over 02- We do substitution a2 = ai + 02 

+00 +00 

4 = ^J(l-2g/a-Z2) f rf^^^-..-lg-2,a, f d5,52e/.+..-l^-,,, 
2crr(— Z2) J J 



SdT{l-2e/a-Z2)ii-'^~^'/'' 
2aT{~-Z2) 



+ 00 

/ daiaj;^""^e"2^"ir(2e/cr + zj, ai^i) , (30) 



where r(A, x) is incomplete gamma function. The value of the last integral can be found in [27j . The final result or 
Ik reads 



h = ^ ^^^ ^,f ^'^ , '^^ r(2e/z2)M-^-/-3-^-/- 2F1 (1, 2e/a, 1 - Z2, 2/3) 
zcr t (t — ^2J 



(31) 



Inserting ([25]) and ([311) into ([Ml) we infer 



'- li^l^^fe^^-'-'"-*-/'' / ^.^.(.^./.-..^/3. 



r(l - 2e/a - Z2) 



T{2TTYT{-e/a-lf'^ 7 (2712)2^ ^v ' / ' - / / .^^ 

— ioo 

r(i + zi) r(-zi + z2) r(£/2 + CTZ2/2) 



v( (^^ M0\v( ( M \mv( /o . ,^rT{-s/2^ az,/2)T{-a{-zi + Z2)/2 - e/2 + a/2) . (32) 

1 (cr(l + zi)/2) 1 {a{-zi + Z2)/2) 1 [a/2 - e - az2/2) 

First we sum over all poles of r(— e/2 — azi/2) and then over pole at Z2 = 0. The result reads 

/2 ^ ^M-^^/^3-/^r(2e/.) 2^.(1, 2s/M,2/3)^?i=^r(./2) 

^ (-1)" r{l-e/a + 2n/a) r{e/a-2n/a) ^, _ 

Z. „, r(.(l-e/a + 2n/.)/2)r(a(e/a-2r./a)/2)'^ e + na/2). (33J 

We will look the final result in the form 

h = M"'"/"e-^^(c_2£"' + c_ie-i) (34) 

To obtain the divergent part of I2 it is convenient to use MATHEMATICA. The result for coefficients c^2 and C-i are 
given by the formula ([TI 



APPENDIX C 

Here we present the derivation of formula (|lip . Expanding scaling function ln(Z) at two- loop order and jl4| one 
can infer that 

~ X! {Kj^'ki + Kk^ji + \i\]k)ninjnknill + - ^ XfjUiUjini + Uj - 2)lf 

l<i<j<k<l<p l<i<j<p 

+ 2_^ i^ij^ik + XijXjk + Xik^jk)ninjnkli ~ 2 2_^ {Xij\k + \jXjk + \kXjk)ninjnkl2 

l<i<j<k<p l<i<j<k<p 

- 22 [KjKknfnjnk + XijXjkntJi^nk + \ik\jkninjn\)ll 

l<i<j<k<p 

- Yl >^ijninj{ni + nj-2)l2 + - ^ X^.^niUjlf . (35) 

By the definition 7 = /i — ^ ' . After differentiation we use the formula \ij(i~^''^ = gjnj + Ag^^Je, which one can 
infer from ([T^. and the integral expansions (|10p and 

n - 1^ + S7^[l^^(4.) + ^(0)(./2)] + 0(eO) (36) 

to derive the following result 

l<z<j<P 1^^<J^P l<i<j<p 

+ ^ n.tijnkigmjgmk + gmjgRjk + gmkgBjk) ( ^ p( /2)2 ln(3/4) 
I] 9%jn,n, {n, + n, - 2) (^^ V{ul2f ^"^^^1^^) ' 



(37) 



l<i<j<p 



The critical exponent is the value of this expression evaluated at the fixed point gj^ij = e. It easy to see that the resut 
is equivalent to (|lll) . 
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